clc;
clear;

index = 1;

syms x y t nu Sc nu_m
switch index
    case 1
    u1 = x.^2.*(x - 1).^2.*y.*(y - 1).*(2*y -1).*cos(t);
    u2 = -y.^2.*(y - 1).^2.*x.*(x - 1).*(2*x - 1).*cos(t);
    B1 = sin(pi*x).*cos(pi*y).*cos(t);
    B2 = -sin(pi*y).*cos(pi*x).*cos(t);
    p = (2*x-1).*(2*y -1).*cos(t);
end

u1dt = diff(u1,t);
u1dx = diff(u1,x);
u1dy = diff(u1,y);
u2dt = diff(u2,t);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
B1dt = diff(B1,t);
B1dx = diff(B1,x);
B1dy = diff(B1,y);
B2dt = diff(B2,t);
B2dx = diff(B2,x);
B2dy = diff(B2,y);
pdx = diff(p,x);
pdy = diff(p,y);

u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
B1dydy = diff(B1dy,y);
B1dydx = diff(B1dy,x);
B2dxdx = diff(B2dx,x);
B2dxdy = diff(B2dx,y);
u1B2dx = diff(u1*B2,x);
u1B2dy = diff(u1*B2,y);
u2B1dx = diff(u2*B1,x);
u2B1dy = diff(u2*B1,y);

f11 = u1dt - nu*(u1dxdx+u1dydy) + (u1*u1dx+u2*u1dy) + pdx + Sc*B2*(B2dx-B1dy);
f12 = u2dt - nu*(u2dxdx+u2dydy) + (u1*u2dx+u2*u2dy) + pdy + Sc*B1*(B1dy-B2dx);
f21 = B1dt + nu_m*(B2dxdy-B1dydy) - (u1B2dy-u2B1dy);
f22 = B2dt + nu_m*(B1dydx-B2dxdx) - (u2B1dx-u1B2dx);

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y,t})
matlabFunction(simplify(u1dy), "Vars", {x,y,t})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y,t})
matlabFunction(simplify(u2dy), "Vars", {x,y,t})
disp("===B1导数===");
matlabFunction(simplify(B1dx), "Vars", {x,y,t})
matlabFunction(simplify(B1dy), "Vars", {x,y,t})
disp("===B2导数===");
matlabFunction(simplify(B2dx), "Vars", {x,y,t})
matlabFunction(simplify(B2dy), "Vars", {x,y,t})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y,t})
matlabFunction(simplify(pdy), "Vars", {x,y,t})
disp("===非定常MHD方程===")
matlabFunction(simplify(f11), "Vars", {x y t nu Sc nu_m})
matlabFunction(simplify(f12), "Vars", {x y t nu Sc nu_m})
matlabFunction(simplify(f21), "Vars", {x y t nu Sc nu_m})
matlabFunction(simplify(f22), "Vars", {x y t nu Sc nu_m})